16 



SDBDDSD 

PARIS RESEARCH LABORATORY 



Contribution 
a la Resolution Numerique 
des Equations de Laplace 
et de la Chaleur 



February 1 992 



Jean Vuillemin 



16 



Contribution 
a la Resolution Numerique 
des Equations de Laplace 
et de la Chaleur 



Jean Vuillemin 



February 1 992 



Publication Notes 



The author can be contacted at the following addresses: 

Jean Vuillemin 

Digital Equipment Corporation 
Paris Research Laboratory 
85, avenue Victor Hugo 
92563 Rueil-Malmaison Cedex 
France 

jv@prl . dec . com 

A video tape exists, demonstrating the described design. Information can be obtained from 
the librarian at PRL. ( librariangprl . dec . com) 



© Digital Equipment Corporation 1992 



This work may not be copied or reproduced in whole or in part for any commercial purpose. Permission 
to copy in whole or in part without payment of fee is granted for non-profit educational and research 
purposes provided that all such whole or partial copies include the following: a notice that such copying 
is by permission of the Paris Research Laboratory of Digital Equipment Centre Technique Europe, in 
Rueil-Malmaison, France; an acknowledgement of the authors and individual contributors to the work; 
and all applicable portions of the copyright notice. Copying, reproducing, or republishing for any other 
purpose shall require a license with payment of fee to the Paris Research Laboratory. All rights reserved. 

ii 



Resume 



Nous adaptons les techniques classiques de resolution par differences finies des equations de la 
chaleur et de Laplace, au traitement par des operateurs materiels specialises. Le parallelisme du 
calcul est obtenu par un pipe-line regulier, dont la profondeur arbitraire s'adapte directement 
au volume de materiel disponible. 

Une implantation de la methode, sur technologie PAM (memoire active programmable 
[BRV 89], [BRV 92]), utilise un pipe-line de 128 operateurs 24 bits, avec une frequence 
d'horloge de 20 MHz. Ceci permet de calculer 5 G 1 operations utiles par seconde, addition ou 
decalage en virgule fixe. Les resultats obtenus en virgule fixe etant, pour ce type de probleme, 
egaux a ceux obtenus en virgule flottante, cette realisation depasse, en valeur absolue, les 
puissances de calculs precedemment publiees ([McB 88] et [McB& al 91]), tous types de 
materiels confondus. Un ordinateur sequentiel classique devra, quant a lui, executer 25 G 
instructions par seconde pour reproduire le meme calcul. Si enfin on compare, en francs par 
operation par seconde, la resolution des memes equations sur super-ordinateur scalaire et 
ordinateur massivement parallele, la technique proposee permet une economie de deux ordres 
de grandeur sur le prix de ce calcul. 



Abstract 



We adapt the classical finite difference method to compute solutions of the Heat and Laplace 
equations with help from special purpose hardware. Parallelism is achieved through a pipe-line, 
whose depth can be adapted to available silicon area. 

An implementation of the method on PAM (Programmable Active Memory) technology runs 
at 20 MHz, with a pipe depth of 128 operators. This lets us compute 5G additions and shifts per 
second, with a 24 bits fixed-point data format. Since it is easy to show that fixed-point gives the 
same results as floating-point for these two problems, this exceeds computing performances 
previously reported ([McB 88] and [McB& al 91]) for super-computers. 

A serial computer will need to execute 25 G instructions per second, to reproduce the same 
computation. The price of solving the Heat and Laplace equations, in $ per operation per 
second, is two order of magnitude lower with the proposed PAM implementation, than with 
super-computers. 



1 IT = 10 12 , 1G = 10 9 , 1M = 10 6 , IK = 10 3 . 
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1 Introduction 

L' analyse numerique, au sens general de [DL 88], est le domaine privilegie d' applications 
des super-ordinateurs. Pour de multiples raisons : 

• La simulation numerique du monde physique est l'une des cles des technologies 
modernes; de la chimie moleculaire au nucleaire, de l'espace a l'electronique. 

• Dans beaucoup de ces domaines, la puissance de calcul aujourd'hui disponible est 
derisoire, en regard de l'ampleur des problemes a traiter; de la simulation d'une aile 
a celle de l'avion tout entier; de la recherche d'une nappe de petrole a la cartographie 
souterraine du globe. 

• De par la nature meme des phenomenes physiques sous-jacents, nombre de ces problemes 
offrent des possibilites, sans limite conceptuelle, au parallelisme des calculs; de la 
simulation des galaxies a celle des circuits integres, de l'equation de Schrodinger a 
celles de Maxwell et Laplace. 

• Au vu de l'importance scientifique, fmanciere ou strategique de certaines applications, 
les agents economiques, etats et industrie, consacrent des sommes considerables a la 
mise au point et a l'exploitation des super-ordinateurs. 

1.1 Super-ordinateurs ou circuits specifiques? 

Pour arriver aux performances des super-ordinateurs, du Gops en 1990 au Tops en l'an 
2000, on suit, pour l'instant, deux voies principales : 

1. les semi-conducteurs (ECL ou AsGa) ultra-rapides; 

2. les structures massivement paralleles. 

Elles sont largement antagonistes. Les technologies rapides sont, aujourd'hui, intrinsequement 
chaudes par leur consommation electrique. Les machines resultantes sont done volumineuses, 
et a faible parallelisme. Les structures massivement paralleles utilisent des technologies 
beaucoup plus lentes. Elles doivent de plus consacrer, en l'etat courant de l'art, une partie 
importante de leurs ressources a la gestion meme du parallelisme. C'est pourquoi il leur faut 
beaucoup de processeurs, pour simplement egaler les performances des machines sequentielles 
les plus rapides. Dans les deux cas, les super-ordinateurs resultants sont d'un usage complexe, 
et coutent tres cher. 

II existe pourtant une troisieme voie, rendue possible par l'avenement des circuits a 
haute densite d' integration (VLSI). Ceux-ci permettent, en theorie, de realiser des materiels 
specifiques, dont la structure physique est directement calquee sur la nature algorithmique de 
chaque application : 
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• Avantages : dans un volume et pour un prix comparable a ceux d'un micro-ordinateur, 
on dispose d'une puissance de calcul comparable a celle d'un super-ordinateur, voire 
nettement superieure dans bien des cas. 

• Inconvenients : une telle machine n'est capable d'executer qu'une classe tres particuliere 
d'applications, voire une seule. Le cout et le temps de realisation du prototype initial 
sont eleves, sans commune mesure avec celui de la premiere mise au point d'un logiciel 
pour super-ordinateur. 

Ces inconvenients, combines aux difncultes que rencontre l'analyste numericien devant 
formuler ses problemes en termes utilisables par un electronicien, cantonnent, pour le moment, 
notre troisieme voie aux etudes academiques, sans retour experimental. Ceci, en depit de son 
impact economique, potentiellement considerable. 

1 .2 Circuits programmables 

La technologie PAM (memoires actives programmables) change les donnees du probleme. 
II s'agit de composants dont on peut changer dynamiquement, plusieurs centaines de fois 
par seconde, la configuration : structure de connexion et tables de verite des fonctions 
booleennes. On consultera les references [GK 89], [X 87] et [K 89] pour des implantations 
et des applications de cette technologie. Nous decrivons les realisations et applications de 
DecPeRLe 0 , une PAM de 50K portes logiques en [BRV 89], et DecPeRLei, une PAM de 
250K portes en [BRV 91]. Ces etudes montrent que : 

1. Tout circuit digital synchrone est realisable, sur une PAM de faille sufflsante. Une 
PAM de faille infinie est au calcul parallele ce qu'une machine de Turing est au calcul 
sequentiel. 

2. Si on compare, a technologie egale, les performances de deux realisations du meme 
algorithme parallele, par un circuit specifique et par PAM : 

• La periode d'horloge est comparable, car invariablement determinee par des 
considerations externes (vitesse des memoires ou du Bus); 

• Le circuit specifique est plus petit, d'environ deux ordres de grandeur, que la 
realisation PAM, soit une carte 100 cm 2 ) pour un composant 1 cm 2 ); 

• Le temps de mise au point d'une application PAM est comparable a celui de 
l'ecriture d'un logiciel tres optimise, quelques semaines, en regard des mois 
necessaires a la realisation d'un circuit specifique; 

• Le cout materiel initial de la solution PAM est nul 2 , celui d'un circuit specifique 
est de quelques centaines de milliers de francs. Le circuit specifique ne devient 
rentable qu'a partir de plusieurs centaines de systemes realises. 

2 a condition, bien entendu, de disposer de cette technologie. 
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1.3 Contribution 

Munis de cette technologie PAM, nous tentons, depuis trois ans, d'en evaluer l'impact 
scientifique, en particulier pour l'analyse numerique. Sous l'impulsion de J.P. Quadrat, nous 
avons choisi de resoudre l'equation de Laplace, en deux dimensions : 

8 2 F | 8 2 F _ 
Sx 2 Sy 2 

C'est, pour J.L. Lions, le probleme fondamental de l'analyse numerique. Un premier essai 
de 40K portes, pour 100 Mops, realise sur DecPeRLeo par F. Rocheteau, nous apporte les 
enseignements suivants : 

• Le schema aux differences finies, en deux dimensions : 

4F n+l (x,y) = F n (x + 1,2/) + F n (x -1,2/) + F n (x,y + 1) + F n (x,y - 1) 

est a abandonner, en faveur du schema, en une dimension, du Paragraphe 4.3. Ce 
dernier divise la bande passante des memoires par deux. Sa topologie lineaire reduit 
notablement la surface de l'operateur atomique, en augmentant d'autant la profondeur 
du pipe-line. Son utilisation, par la methode des directions alternees, permet de resoudre 
a la fois l'equation de la chaleur et celle de Laplace, en dimension arbitraire, et non 
seulement deux. 

• L'utilisation d'une arithmetique entiere, plutot que flottante, est la cle des performances 
de ce type de materiel; la precision initialement choisie de 16 bits n'est cependant pas 
suffisante, pour les tailles de grilles compatibles avec la vitesse du traitement; il faut 
passer a 24 bits de precision ou plus. 

• Le prix a payer, pour simplifier l'operateur atomique, est un pas de temps tres petit, 
fixe par l'equation (4). Ceci est bon pour l'equation de la chaleur, mais trop lent pour 
celle de Laplace. On accelere la convergence, dans ce cas, par l'emploi des techniques 
multi-grilles du Paragraphe 6.3. 

Nous justifions ici les decisions d'implantation, et decrivons les resultats obtenus. Les 
Paragraphes 4.4.1 (plomberie numerique), et 6.3 (multi-dimensions), sont les seuls oil nous 
nous eloignons de methodes bien connues, en analyse numerique; pour tous les autres usages 
de techniques classiques (et largement elementaires), nous renvoyons le lecteur aux ouvrages 
de base, principalement [DL 88]. 
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2 Equation de la chaleur 

On se propose de resoudre numeriquement l'equation de diffusion de la chaleur : 

6T 

valable dans un domaine CI C Ck interieur au cube Ck = [0, l[ k de l'espace R k de dimension 
k > 0 (typiquement k = 2, 3). Dans l'equation (1), le symbole A represente l'operateur de 
Laplace (on dit aussi laplacien) : 

AT = — + ••• + — 

8x\ bx\ 

Les temperatures initiales, au temps t = 0, sont donnees, en tout point 
(x\ , • • • , Xk) e Ck du cube, par la condition de Cauchy : 

T(0,x u ---,x k ) = T 0 (x u ---,xk). (2) 

On suppose aussi connues, et fixes Vf > 0, les temperatures en tout point (xi , • • • , Xk) e T du 
complementaire de Q. dans Ck = Q U T, par la condition de Dirichlet : 

T(i,a:i, • • -,x k ) = T 0 (x u - ■ ■ ,x k ). (3) 

On conviendra, sans perdre en generalite pour la suite, que T est compose d'un ensemble fini 
de sources ponctuelles de chaleur (representees mathematiquement par des masses de Dirac, 
au sens des distributions). On etend enfm les solutions de (1,2,3) a tout point de l'espace R k 
par periodicite, soit, pour 1 < j < k : 

T[xj + 1] = T[xj]. 

Dautray, Lions [DL 88] et Feynman [FLS 63] decrivent quelques unes des multiples appli- 
cations du calcul de (1,2,3) en thermodynamique, neutronique, cinetique chimique, hydrody- 
namique, . . . 

Remarque : En l'absence de la condition (3), soit r = 0, £1 = Ck, les solutions de (1,2) sont 
telles que : 

S f 

— T(t,x l ,---,Xk)dx l ---dx k -0, 
bt J Ck 

pour t > 0, ce qui exprime la conservation de la quantite initiale de chaleur. 



3 Resolution par differences finies 

On se fixe un pas d'espace Aa; = Axi - ■ ■ ■ = Axk, identique dans chaque dimension, et un 
pas de temps At tel que : 

k 

At = -Ax 2 . (4) 
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On calcule les valeurs approchees G n {i\ , • • • ,ik) des solutions T(t n , x^ , ■ ■ ■ , Xi k ) de (1,2,3), 
prises en t n = nAt, X{ = iAx, par un schema (explicite) aux differences finies, que Ton peut 
ecrire (methode des directions alternees) : 

G n +i(i\, ■ ■ -,ik) = Hi ■ ■ -UkGndi, ■ ■ -,ik) (5) 

La relation (4) permet de donner a chacun des operateurs Tij, pour 1 < j < k, la forme 
particulierement simple : 

HjG n (ii,---,ik)- si (x { , ■ ■ -,x k ) e F alors G n (i\, ■ ■ ■ , i k ) 

sinon j(G[ij + 1] + 2G + G[ij - 1]). 1 ' 

La formule (5) ne depend pas du choix du syteme d'axes, puisque les operateurs (6) commutent 
tous, deux a deux : 

TijTii — 7i{H.j. 

Cette propriete, dont nous ferons bon usage, traduit le caractere isotrope de la diffusion de la 
chaleur, dans le monde physique. 

Remarque : Quand r = 0, le schema (5) conserve la quantite initiale de chaleur, pour tout 
n > 0 : 

^2 G n (x { , ■ ■ -,x k ) = ^2 Gq(x { ,- ■ ■ ,x k ). (7) 



4 Structure materielle 



Dans ce qui suit, nous decrivons une structure materielle permettant de tirer le meilleur parti 
des technologies de circuits VLSI, dans le calcul des solutions de l'equation de la chaleur. Ann 
d'alleger les notations, on se place en dimension k = 2; la generalisation a k > 2 est traitee au 
Paragraphe 4.5. On se propose done de resoudre numeriquement le probleme : 



T(0,x,y)=T 0 (x,y) 



st 

ST. 
St 



T(t,x,y) = 0 



S 2 T + S 2 T 
Sx 1 Sy 2 



T(t, x,y) = T(t, x + l,y) = T(t, x,y+l) 

On utilise pour cela le schema aux differences finies : 

G n +i(i,j) = H x H y G n {i,j), 

n x G(i,j) = kG(i + l,j) + 2G + G(i- 

n y G(i,j) = \(G(i,j + \) + 2G + G{i,j 

G(i,j) = n x G(i,j) = H y G(i,j), 



pour 0 < x , y < 1 , 
pour x, y g r, t > 0, 
pour x, y e r, t > 0, 
pour t > 0. 



- D), 



pour i, j e T, 
pour i, j e T, 
pour i, j g r. 



(8) 



(9) 



4.1 Temperatures discretes 



Soient / = inf Tq et S = sup Tq les valeurs extremes des temperatures initiales dans le 
carre C*2. Par le principe du maximum, la solution de (8) satisfait 

I <T(t,x,y)< S, 
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pour tout t > 0 etO < x,y < \; ceci est evident si Ton se reporte a l'interpretation physique 
du probleme. On represente alors les temperatures discretes G n (i, j) par des entiers naturels, 
avec la correspondance : 

5-7 

T(nAt, iAxJAy) = I + 2b _ { G n (i, j). 

Le nombre b de bits dans l'ecriture binaire de ces entiers est tel que : 

2~ b < sup \T Q (iAx + e, jAy + e) - T Q (iAx, jAy)\/T 0 (iAx, jAy). 

De cette faqon, les erreurs dues a l'arrondi numerique sont inferieures a celles inherentes a 
la discretisation des temperatures initiales. Avec ce choix, les decalages lies aux exposants, 
dans une representation des temperatures par des nombres flottants, sont factorises une fois 
pour toutes, a la presentation initiale du probleme. Tous les calculs ulterieurs sont faits sur la 
mantisse entiere, ceci en vue de simpliner, sans perte de precision numerique, la realisation 
materielle de l'operateur H, qui suit. 



4.2 Representation du probleme 



La grille de travail est de taille m x m, l'entier m = 2 a etant une puissance de deux. Les 
pas correspondants d'espace et de temps sont donnes par : 

Ax = Ay = 2~ a , At = 2~ (2a+l) . 

A l'instant t n = nAt, l'etat du systeme est stocke dans une memoire composee de m x m 
mots contigiis, de b + 1 bits chacun. Aux b bits representant la temperature discrete, nous en 
ajoutons un, qui vaut 1 si le point correspondant de la grille est une condition limite 
(iAx,jAy) e T, sinon 0. En convenant de le placer en tete des poids faibles de l'ecriture 
binaire, on reduit le test d'appartenance a T, a un test de parite, et on trouve, a l'adresse i + mj 
du tableau M, la quantite : 



M(i + mj) = 



l + 2G n (i,j) si(i,j)eT, 
2G n (i,j) si(i,j)er. 



La memoire M est dotee d'un controleur, generant des suites d'adresses, qui permettent de 
parcourir son contenu, de deux manieres : 



1. suivant l'axe des x, dans un balayage horizontal, correspondant a l'ordre lexicographique 
(h j) <x (i', j') si i < i', ou i = i' et j < j'; 

2. suivant l'axe des y, dans un balayage vertical, correspondant a l'ordre lexicographique 
dual (i, j) < y (i', j') si j < j', ou j = j' et i < i'. 



Le controleur memoire est compose de deux compteurs binaires sur 2a bits, et d'un multi- 
plexeur : 
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• Lecompteur horizontal C^genereles 2a bitsd'adresse dans l'ordre&o ' ' ' b a -ib a • • • &2a-i; 

• Le compteur vertical C y inverse les a bits de poids fort, avec ceux de poids faible, dans 
l'ordre b a ■ ■■b 2a -ih) ■ • -6 a -i; 

• Le multiplexeur connecte les adresses de la memoire, a l'un ou l'autre des deux 
compteurs, suivant la nature de la passe. 

Remarque : Avec un tel controleur, la structure topologique de l'espace de travail est 
naturellement celle d'un tore a deux dimensions. En remplacant la condition de periodicite : 

T(t,x,y)=T(t,x + l,y) = T(t,x,y+l), 

par une condition de Neuman sur les bords du carre : 
8 8 

—T(t,x,b)= —T(t,b,y), pour b = 0 ou b = 1, 
8x 8y 

on retrouve la topologie usuelle du plan euclidien. Les modifications correspondantes, dans le 
controle des compteurs C x et C y , sont les suivantes : 

1. lors d'un balayage horizontal, la valeur de C x doit etre tenue constante (en fixant son 
increment a zero), pendant un nombre de cycles egal a la profondeur p du pipe-line des 
donnees (cf. Paragraphe 4.3), avant (C x = - 1 (mod m)), et apres (C x = 0 (mod m)) 
chaque retour de ligne; 

2. lors d'un balayage vertical, la valeur de C y doit etre tenue constante, pendant p cycles, 
avant (C y = - 1 (mod m)), et apres (C y = 0 (mod m)) chaque retour de colonne. 

4.3 Traitement a la chaTne 

Pour resoudre les equations (8), nous nous dotons de : 

• Une memoire M. s , munie de son generateur d' adresses, servant de source des donnees. 
Elle represente l'etat du systeme au temps t n . 

• Une memoire Md, servant de destination aux resultats du calcul. Elle est identique 
a M. s , et contiendra, en fin de passe, l'etat du systeme au temps t n+ \. La transition 
TjHn + 1 inverse les roles de M s et Md- On peut, sans obstacle theorique, representer 
M s et Md dans la meme memoire physique; ceci divise cependant la bande passante 
du traitement, par un facteur deux, puisqu'il faut lire et ecrire une memoire a chaque pas 
de calcul. 

• Un operateur 7i, dont on trouvera la description au Paragraphe 4.4; son objet est de 
transformer les donnees provenant de M s en donnees a destination de Md- 
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• Un controle externe, qui etablit le sens du transfert M s ^> Md, et la nature de la passe : 
horizontale ou verticale, avec ou sans multi-dimension (cf. Paragraphe 6.3). 

Voici, sous ces hypotheses, le schema d'une iteration (9), composee d'une passe horizontale 
suivie d'une passe verticale a travers tout l'espace : 



M 


i— > 


H 


i— > 

X 


M 


i— > 


H 


i— > 

y 


M 



II faut repeter n fois cette operation, arm d'obtenir (7i x 7i y ) n Go 
diffusion, nous avons : 

(7i x 7iy) p G r = 7i x 7iyG r = G r+P . 

Ainsi, le calcul de p passes : 



G n . Par isotropic de la 



M 


i— > 


H 


h- 

X 


■* M 


1— » 


H 


h- 

y 


-» M 


\—> 


i— > 


M 


i— > 


H 


i— » 

X 


M 


i— » 


H 


y 


M 



donne le meme resultat qu'une seule passe de p traitements Ti a la chaine : 

v v 



M 


i— > 


H 




H 


i— > 


M 




H 


y 


H 


i— > 


M 



Ceci justifie d'utiliser le schema : 

G n+P (i,j) = H p x H p y G n (i,j), 



(10) 



qui cumule p passes en une, avec exactement la meme demande de bande passante memoire 
que le schema (9), pour une seule passe. En choisissant la plus grande valeur de p realisable 
dans une technologie donnee, on adapte le calcul de chaque iteration au volume de logique 
disponible, en deroulant la boucle interne du calcul, autant que faire se peut. La structure 
d'ensemble du systeme est done : 



H 



H 



M d 
C d 



(11) 



Remarque : Le pipe-line des donnees introduit un retard de p cycles entre les adresses des 
sources M s , et celles des destinations Md', on compense en retardant de p cycles la mise en 
route des compteursdeA^d,afmdemaintenirinvariantes les relations : C£ = C*+p, Cy = C*+P- 



4.4 Operateur atomique 

L' ensemble du systeme (11) est synchronise par une horloge digitale r, de periode At. 
Comme une etape du schema (10) correspond a deux parcours complets des memoires (en x 
et en y), on a la relation suivante entre le pas de calcul At et le pas de temps At : 

m At - pAt. 
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Dans ces schemas, les carres sont des registres synchrones 1 bit (flip-flop), les petits losanges 
des multiplexeurs, et les grands des additionneurs binaires complets. La temperature d' entree 
est T = £ ; > 0 (T; + Tj)2\ et celle de sortie: 

HT = Y^(HTi + HT!)2\ 

i>0 

Figure 1: Schemas de 7i = 7i~7i + , sur 4 bits. 
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A A' B B' A A' B B' A 

i-1 i-1 i-1 i-1 i i i i i- 



I I J 



S S' S S' S S' 

i-1 i i i+i i + 1 i+2 

Ici, les losanges representent des additionneurs binaires complets. Les temperatures d'entree 
sont T a = Ei>o(^i + 4)2*, T b = E*>o( 5 i + 5 ;)2\ et celle de sortie T s = T a + T b = 
Ei>o(Si + Sl)2\ 

Figure 2: Addition sans retenues, sur 2 bits. 

Dans une passe horizontale (resp. verticale), les entrees des operateurs 7i arrivent dans l'ordre 
< x (resp. < y ); chaque operateur 7i calcule un pas de diffusion en une dimension, suivant 
l'axe determine par l'ordre d'arrivee de ses entrees. On decompose 7i = 7i + 7i~ = 7i~7i + en 
produit d'un operateur 7i + de diffusion en avant, par un operateur 7i~ de diffusion en arriere. 



H + comprend un registre ^pj et un circuit combinatoire A + 



, connectes comme suit : 



G(v + 1) ^ [7] ^ G(v) 
G(v + 1) 



7i + G(v). 



H comprend un registre et un circuit combinatoire _A_, connectes comme suit : 

► n~G(v). 



G(v) ^ 0 h* G(v - 1) 
G(v) 



A~ 



Si E = (E1E2 • • • Ek • • *) represente la suite des valeurs d'entrees d'un registre 

E ^ g S, 

la suite des sorties S = (OE1E2 ■ ■ ■ Ek-i ■ ■ •) correspondantes est la meme, a un retard d'une 
periode At pres. Les circuits A~ et A + sont identiques, apres permutation des entrees : 



A(v 0 ,vi) = A + (v 0 ,vi) = A (v u v 0 ), 

'•• ' " 1 si v 0 € r, 



A(v 0 ,vi)= jiVo + n), 
A(v 0 ,vi) = v 0 , 



si g r. 
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4.4. 1 Plomberie numerique 

Afin de fixer la structure de l'operateur A, il convient de regler la question d'arrondi liee a 
la division par deux ci-dessus. Le choix naturel : 

A(v 0 , vi) = [_ — - — J , si v 0 e r, 

amene des pertes numeriques de chaleur en violation de 1' invariant (7). Cette question, 
qui n'est pas abordee dans notre bibliographie, a pourtant une incidence significative sur les 
resultats numeriques, des que la taille m 2 de la grille devient grande, en regard du nombre b 
de bits representant les temperatures. La solution la plus simple a ce probleme est d'interdir le 
partage du bit de poids faible des temperatures (de poids 2 puisque le bit de parite represente 
l'appartenance a F). Ceci est realise par la formule arithmetique : 

f A(v 0 ,vi) = v 0 si v 0 -|- 2= 1, 

I A(v 0 , vi) = 2(v Q + 4 + v\ + 4 + (v 0 | 4) + 2), sinon. 

Ici, n-\ - d designe le reste de la division entiere de n par d, et n d le quotient. La formule 
(12) restaure l'invariant (7), de la quantite initiale de chaleur. 

4.4.2 Additions sans retenues 

Le pas de calcul At est determine par le composant le plus lent (chemin critique) de 
l'ensemble du systeme synchrone (1 1). Notre objectif est de regler At sur la vitesse nominale 
des memoires. Nous montrons en [V 91] comment realiser des compteurs binaires capables 
d'operer a la vitesse intrinseque (toggle rate) de la technologie de realisation, pour un nombre 
de bits arbitraire. Ceci resoud la question des controleurs d'adresse. 

Pour regler celle des additionneurs A, on utilise la technique classique (carry save adder) 
des representations binaires redondantes. A l'interieur du chemin des donnees, chaque 
temperature T est representee par la somme T = R + S de deux nombres binaires R et S. 
L' addition de deux temperatures redondantes T + T' = R + R' + S + S' se calcule alors bit a 
bit, sans propagation de retenue, comme dans les sche mas d e la Figure 2, qui sont realises 



exclusivement a partir d' additionneurs binaires complets | abc 
Le circuit combinatoire 



a 
b 
c 



abc 



calcule l'unique solution binaire a, b, c,r,ss {0, 1} de l'equation : 

a + b + c-s + 2r. 



Avec cette technique, le delai critique de l'operateur H est double de celui de | abc 
independamment du nombre de bits b; on le rend facilement inferieur a celui des memoires. 
L' inconvenient de la methode est de doubler la surface necessaire a la realisation de H, 



Research Report No. 16 



February 1992 



12 



Jean Vuillemin 



en comparaison d'une propagation combinatoire des retenues sur b bits (voir la remarque 
ci-dessous). 

En fin du chemin des donnees, apres les p operateurs 7i, on doit ajouter un unique 
additionneur combinatoire (voir par exemple [GV 82] pour une structure de delai 2log2b), 
qui convertit la representation redondante R + S en un entier binaire unique T = R + S, 
pour ecriture dans la memoire M.d- On rendra ce circuit plus rapide que les memoires, en 
coupant son chemin critique par un nombre suffisant de registres binaires (au plus log2b). La 
surface totale de cet unique element est negligeable devant celle des p operateurs 7i, pour p 
suffisamment grand; le nombre r de registres, mis en serie dans cet additionneur, augmente 
d'autant la longueur totale du pipe-line, qui devient p + r. 

Remarque : En pratique, il est inutile de pousser la representation redondante jusqu'au 
niveau du bit. Decoupons, en effet, b - cden d blocs de c bits; la propagation des retenues est 
combinatoire, sur c bits, a l'interieur de chaque bloc; d bits redondants suffisent alors a casser 
les chames de retenues, en choisissant pour c la valeur maximale compatible avec le debit des 
memoires. Cette remarque reduit la surface totale d'un tel operateur par un facteur proche de 
deux, sans perte en vitesse pour l'ensemble du systeme. Les schemas de la Figure 1 montrent 
une realisation, sur c = 4 bits; ils illustrent aussi la realisation de la formule (12). 

4.5 Autres dimensions 

Pour appliquer notre systeme a la resolution de l'equation de la chaleur (1,2,3) en k > 2 
dimensions, il suffit d'adapter les controleurs des memoires; le chemin des donnees (forme de 
p operateurs 7i) reste inchange. Chaque controleur comprend alors k compteurs sur a = ka! 
bits. Ainsi, en k - 3 dimensions, on arrange les 3a bits de chaque compteur dans l'ordre : 

• b 0 ■ ■ - b a i_ib a i ■ ■ -&2a'-i&2a' ' ' 'Ha 1 -! pourC^; 

• &a'--- b 2a '- 1 ha' ■ ■ ■ ha'- 1 &0 • • • K'- 1 P our Cy \ 

• ha' ■ ■ ■ ha'- 1 h ■ ■ ■ b a >- 1 b a i ■ ■ ■ b 2a '~ i pour C z . 

On obtiendra, bien entendu, exactement le meme resultat en echangeant arbitrairement les a' 
bits de poids fort de l'un de ces compteurs, avec les a' bits du milieu. La generalisation a 
k > 3 est directe, bien que peu realiste, vu les tailles presentes des memoires; la methode des 
elements finis prend ici le pas sur celle des differences finies. 
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5 Equation de Laplace 

Le cas statique de l'equation de la chaleur, dite equation de Laplace (aussi de d'Alembert), 
est particulierement interessant, par le nombre de lois de la physique qui se reduisent a 
l'equation : 

AU = 0, (13) 

sujette a la condition de Dirichlet : 

U(x x , ■ ■ -,x k ) = U 0 (xi, ■ ■ -,x k ), (14) 

valable en tout point (x\, ■ ■ - ,Xk) e T d'un ensemble discret F. Nous resolvons (13,14) 
en prenant le cas limite, pour t i-> °o, de l'equation de la chaleur (1,2,3). Soit, en effet, 
T(t, xi , • • • , Xk) une solution a l'equation (1) de la chaleur, obtenue a partir d'un etat initial (2) 
arbitraire, compatible avec la condition de Dirichlet (3), que Ton prend identique a (14). La 
fonction 

U(x u - ■■,x k ) = T(oo,X { , ■ ■ -,x k ) 

est solution des equations (13,14); cette solution est unique, independamment de l'etat initial, 
pourvu que T ^ 0 soit non vide. 

Dautray, Lions [DL 88] et Feynman [FLS 63] decrivent quelques applications du calcul de 
(13,14) en electrostatique, electromagnetisme, optique, elasticite, mecanique, hydrostatique, 
etc 

6 Vitesse de convergence 

L'iteration (5) Hi ■ ■ - TikGt =>■ G t +\ converge, pour t i— > oo, vers une solution F - G«, 
independante de etat initial Go, pourvu qu'il soit compatible avec la condition de Dirichlet 
(14), supposee non vide. Le choix de l'etat initial affecte cependant la vitesse de convergence, 
i.e. le nombre minimal T = T(e) de pas d'iteration (t i— > t + 1) necessaire pour obtenir une 
solution Gt de (13,14), approchee a e > 0 pres : 

\AF\ « \G T+ i - G T \ < e. 

En general, le processus devient numeriquement stable, a e = Aa; pres, apres un nombre de 
pas d'iteration (5) de l'ordre de m = Ax~ l . Le nombre total d'operations necessaires a la 
resolution du probleme (13,14) de Laplace, en prenant la limite de revolution des equations 
(1,2,3) de la chaleur, est done proportionnel a m k+l . Les techniques qui suivent permettent de 
reduire ce temps de convergence, par un choix approprie de l'etat initial. 

On se limite ici, pour simphner les notations, au probleme de Laplace en deux dimensions, 
sur une grille de taille m 2 . La complexite du processus iteratif est 0(m 3 ), dans ce cas. 
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6.1 Multi-grilles 

L'idee generate est de choisir, pour etat initial sur une grille fine m x m, une injection 
de la solution au meme probleme de Laplace, calculee recursivement sur une grille grossiere 
m' x m', avec m' < m. La description la plus simple de cet algorithme vient en prenant 
m = 2m! : 

1. Les m 2 donnees initiales Go(i,j) sont projetees sur la grille m' 2 = z j : , suivant la 
formule : 

n ,,., */\ Go(i,j) + G 0 (i + l,j) + Go(i,j + l) + G 0 (i + l,j + l) 
Gq(i , J ) = , 

avec i = 2i' et j = 2j'. On decide de l'appartenance (i', j') e V, aux points de Dirichlet 
de la grille grossiere, par : (2i', 2j') e Tou (2i' + l,2j') e Tou (2i', 2j' + 1) e T, ou 
enfin(2i'+ 1,2/+ 1) e T. 

2. On resoud recursivement le probleme de Laplace, sur les m' 2 points de la grille grossiere. 
La recursion se termine par un calcul iteratif, quand la taille de grille est suffisamment 
petite : m < mo. 

3. On injecte la solution G'„(i' , j') ainsi obtenue pour la grille m' 2 , sur la grille fine : 

Gi(i, j) = G'Ji -f2j'i 2) pour i, j £ F, 
Gi(i, j) = G Q (i, j) pour i, j e T. 

4. On termine le calcul iterativement, en appliquant le schema (9) jusqu'a la convergence, 
a partir de l'etat initial G\(i, j) calcule a l'etape precedente. 

Une analyse simple montre que le nombre total d'operations |m 2 = m 2 + ^m 2 + j^m 2 + • • • 
de cet algorithme, est proportionnel a la taille m 2 du probleme initial. 

6.2 Multi-echelles 

L algorithme multi-grilles s'adapte mal au materiel presente en Section 4 : 

• Le temps de projection (etape 1) et d'injection (etape 3) est egal a celui m 2 Ar d'un 
parcours complet de la memoire; ces operations necessitent, de plus, un materiel 
specifique. 

• La resolution recursive (etape 2) necessite log2m compteurs d'adresse, nombre qui croit 
arbitrairement avec la taille m de la grille. 

L' algorithme multi-echelles (parallel superconvergent multigrid de [FMcB 88]) elimine la 
projection et l'injection du multi-grilles, en resolvant recursivement quatre problemes de 
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Laplace, au lieu d'un seul pour le multi-grilles, sur les grilles grossieres (0,0), (0, 1), (1,0) 
et (1, 1), obtenues en considerant (mod 2). Cette methode converge plus vite (voir 

[McB&al 90]) que le multi-grilles. Chaque etape recursive traite un volume de donnees 
egal a celui m 2 du probleme initial, pour un total proportionnel a m 2 log2m operations. Cet 
algorithme est done plus lent que le precedent, sur des machines sequentielles, ou a faible 
parallelisme. II se montre pourtant plus rapide, sur des machines hautement paralleles, qui 
tirent parti des avantages que le multi-echelles possede sur le multi-grilles : 

• nombre total de passes plus faible, du a sa meilleure convergence; 

• simplicite du controle, en particulier lors des transitions entre passes; 

• volume constant des calculs pour chaque passe, qui permet d'exploiter un parallelisme 
arbitraire. 

II ne reste plus, pour adapter l'algorithme multi-echelles au materiel decrit en Section 4, qu'a 
en eliminer la recursivite. 

6.3 Multi-dimensions 

Si on elimine la recursion, le multi-echelles revient a resoudre, d'abord, un probleme de 
Laplace sur le parallelepipede | x | X 4, en trois dimensions; le resultat de ce calcul est alors 
injecte comme valeur initiale d'une resolution iterative sur le carre m 2 . C'est la structure de 
l'algorithme multi-dimensions qui suit, apres optimisation du choix des dimensions, dans la 
premiere etape. Pour le decrire, remarquons qu'il existe une bijection naturelle entre le carre 
m 2 et le cube en quatre dimensions ^fm A ' . En posant m = 2 2a , done r = ^fm = 2 a , le point 
generique C2 du carre m 2 est mis en correspondance avec celui (u, v, w, x) e C4 du 

cube r 4 par la formule : 

i — u + rv, j — w + rx; 

et reciproquement par : 

u = i -\-r, v = i r, w = j -\- r, x = j + r. 
Le probleme de Laplace, en dimension deux, se resoud alors en deux passes : 

1. On calcule iterativement revolution de l'equation de la chaleur en quatre dimensions, 
pour un nombre d'etapes de l'ordre de *Jm\ le cube r A de depart est en correspondance 
naturelle avec la grille initiale m 2 , pour m = 2 2a , r = *Jm = 2 a . 

2. On calcule iterativement revolution de l'equation de la chaleur en deux dimensions, a 
partir de l'etat calcule a l'etape precedente, par la correspondance naturelle. 

Nous conjecturons que le nombre de passes necessaires a l'etape 2 du multi-dimensions est 
de l'ordre de ^fm, ce qui donnerait a cet algorithme une complexite en 0(vn}^fm). Cette 
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Figure 3: Structure du systeme. 



hypothese est en accord avec les essais experimentaux, ainsi que les analyses que nous avons 
pu mener, dans quelques cas particuliers. 

Remarque : On peut, bien entendu, iterer plusieurs fois la technique multi-dimensions; a 
la limite, on retrouve le multi-echelles. Ceci n'est cependant pas necessaire, vu les tailles de 
grilles realisables, avec les memoires actuelles. 

7 Realisations concretes 

Nous presentons ici quelques caracteristiques des materiels que nous avons realises, au 
moyen de la technologie PAM. Si l'utilisation de composants configurables s'est revelee 
indispensable dans la mise au point des compromis logiciel/materiel auxquels nous sommes 
arrives, les resultats de cette etude sont directement applicables aux technologies, plus 
classiques, de circuits prediffuses (ASIC gate array) et a facon (full custom VLSI). L'utilisation 
de circuits specifiques permettra de reduire, a performance constante, le volume et le prix 
unitaire du materiel necessaire. 

L' ensemble du systeme (hors memoires) correspond a un prediffuse de 180K portes, ou un 
circuit a fagon de 800K transistors, dote de 256 pattes de communication. 

Les temperatures sont representees sur 25 bits, soit 24 bits utiles, et un bit de parite 
pour decider de l'appartenance a Y. Ce choix permet des largeurs m de grilles jusqu'a 
m = 2 24 w 8M, suffisant en pratique. 

7.1 Structure d'ensemble 

Le systeme de la Figure 3, se compose de : 

• Un ordinateur hote (DS5000, station 25 MIPS) qui controle l'ensemble (nature et sens 
des passes). II permet aussi la visualisation en temps reel (de 5 a 20 images 512 x 512 x 8 
par seconde) de l'etat interne du systeme. En mode 3D, les RAM de l'hote servent de 
memoire de donnees, ce qui permet de traiter des cubes de 256 3 echantillons. 
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• L'hote est relie a la PAM par un lien TurboChannel, a haut debit (100M octets par 
seconde). 

• La partie PAM, realisee sur DecPeRLei, comprend : 

- 4M octets de RAM, utilises en mode 2D pour stacker quatre images de 512 2 
echantillons. 

- Un generateur d'adresses, qui comprend les 8 compteurs necessaires au multi- 
echelle, en 2D. 

- Un chemin des donnees, compose de deux pipe-lines de 64 operateurs 7i chacun. 
Chaque element 7i opere sur des entiers 25 bits, organises en 6 blocs de 4 bits, avec 
un bit de redondance par bloc; le bit de parite (appartenance a T) est duplique, par 
raison de symetrie, pour un total de 32 bits par operateur. Cet ensemble correspond 
a 150K portes logiques, soit plus de 80% du total. 

- Un aiguilleur sert a selectionner la source des donnees (note en 3D et PAM en 
2D), et a connecter les deux moities du pipe-line, en parallele, ou en serie. II 
comprend en outre les deux additionneurs, qui assurent la conversion en format 
non redondant, a la fin des pipe-lines. 

7.2 Performances 

Tout le systeme opere a la frequence, 20MHz, des memoires, pour calculer 128 operations 
H, a chaque pas : At = 50ns. 

En optimisant le code machine, pour le jeu d' instructions MIPS, nous arrivons a un total 
de 10 instructions par operation H, compte tenu des conditions de Dirichlet, de la formule 
(12), et de la gestion des adresses. Un ordinateur sequentiel classique devra done executer 25 
G instructions par seconde, pour reproduire le calcul de notre materiel. Notre co-processeur 
PAM, double le volume de son hote (DecPeRLei se loge dans la meme boite que le DS5000), 
et multiplie ses performances par mille, pour ce probleme particulier. 

En ne comptant, maintenant, que l'arithmetique dans (12), on arrive a 5 G operations 
effectives par seconde, addition ou decalage en virgule fixe; avec, en prime, conservation de 
la quantite initiale de chaleur (7), et prise en compte des conditions de Dirichlet. Ceci depasse 
les puissances de calculs mesurees par [McB 88] et [McB&al 90], qui donnent la Connection 
Machine CM-2 a 3.8Gflop, et le CRAY Y-MP a 0.7Gfiop. 

D'apres [BKH 88], la mesure, en $ par Mflop, de ces super-ordinateurs est de 3K$/Mflop 
pour la Connection Machine CM-2, et 10K$/Mflop pour le CRAY Y-MP. En estimant, sur la 
base (discutable) des 20 prototypes construits, le prix de PeRLe\ a 100K$, nous trouvons 
25$/Mop, soit un gain de deux ordres de grandeur sur le prix de la resolution des equations de 
Laplace et de la chaleur. 
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8 Conclusion 

Nous n'avons, bien entendu, aborde qu'un aspect partiel, et tres simple, de l'analyse 
numerique contemporaine. II reste a comprendre, et c'est loin d'etre immediat, comment 
appliquer cette technologie a des methodes plus generales (diffusion avec second membre, 
milieu anisotrope, elements finis), et aux autres equations de la physique (ondes, transport, 
Navier-Stokes, Schrodinger). 

Les resultats obtenus sont cependant positifs : ils tracent un chemin direct menant a des 
solveurs numeriques specialises, dont l'unite de puissance est le Top/sec, mille fois plus rapide 
que les super-ordinateurs actuels. C'est techniquement faisable, en quelques annees, et sur un 
petit budget. 

Les analystes numericiens sauront-ils exploiter des machines dont les possibilites algorith- 
miques sont aussi contraintes que celles que nous proposons ? 
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